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Abstract 

The density functional approach in the Kohn-Sham approximation is widely used to 
study properties of many-electron systems. Due to the nonlinearity of the Kohn-Sham 
equations, the general self-consistence searching method involves iterations with alternate 
solving of the Poisson and Schrodinger equations. One of problems of such an approach 
is that the charge distribution renewed by means of the Schrodinger equation solution 
does not conform to boundary conditions of Poisson equation for Coulomb potential. The 
resulting instability or even divergence of iterations manifests itself most appreciably in 
the case of infinitely extended systems. The published attempts to deal with this problem 
are reduced in fact to abandoning the original iterative method and replacing it with some 
approximate calculation scheme, which is usually semi-empirical and does not permit to 
evaluate the extent of deviation from the exact solution. In this work, we realize the itera- 
tive scheme of solving the Kohn-Sham equations for extended systems with inhomogeneous 
electron gas, which is based on eliminating the long-range character of Coulomb interac- 
tion as the cause of tight coupling between charge distribution and boundary conditions. 
The suggested algorithm is employed to calculate energy spectrum, self-consistent poten- 
tial, and electrostatic capacitance of the semi-infinite degenerate electron gas bounded by 
infinitely high barrier, as well as the work function and surface energy of simple metals in 
the model with homogeneous distribution of positive background. The difference between 
self-consistent Hartree solutions and those taking into account the exchange-correlation 
interaction is analyzed. The comparison with results previously published in the litera- 
ture is carried out. The case study of the metal-semiconductor tunnel contact shows this 
method being applied to an infinitely extended system where the steady-state current can 
flow. 

PACS:71.15.-m, 71.15.Mb, 65.40.gh, 85.30.Mn 

1 Introduction 

In the density functional approach, the set of equations describing the inhomogeneous electron 
gas in the Kohn-Sham approximation should be satisfied by the self-consistent distribution of 
electron density A^(r) and Coulomb potential U{r) Due to the essential nonlinearity of 
these equations, the only general method of self-consistent solution involves iterations using 
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alternately the Poisson equation for potential and the set of Schrodinger equations for single- 
particle wave functions in the effective potential Ucs{r) = U{r) + t/xc(r), where U^c is the 
exchange-correlation potential in the local density approximation. A well-known problem of 
such an approach lies in the necessity of taking into account the boundary conditions of the 
Poisson equation, which impose certain requirements in terms of integral relations on the charge 
distribution in the right-hand side of this equation [2]. From the physical point of view, these 
integral relations specify, for example, the total charge of the electron system, if the boundary 
conditions are imposed on the electric field, or the total dipole moment of the system, if the 
values of potential on the boundaries of an infinite region are fixed. 

However, until the self-consistency is reached, the distribution of electron density found by 
means of the Schrodinger equation at each step of iteration process turns out, as a rule, to be 
incompatible with the Poisson boundary conditions specified. Consequently, as shown in j2] , in 
the case of infinitely extended systems it is either not possible to build the next step solution 
at all, or the iteration process lacks convergence as reported, for example, in [3l E]. In such 
a situation, some authors resort to replacing the process of solving the Poisson equation with 
some kind of variation scheme relative to numerical parameters that specify the functional form 
chosen to approximate the potential or charge distribution. This is discussed in more detail in 
subsections 13.11 and 14.11 

To manage this difficulty there are several artificial techniques suggested, all of them chang- 
ing the parameters of the system in question, for example, the compensating background charge 
density [S]. However, the existence and uniqueness of the solution of self-consistent field equa- 
tions follows from their being obtained as Euler-Lagrange equations for a variational problem 
of minimizing the total energy of non-degenerate ground state (see, for example, [1]). In such a 
case the single-particle orbital wave functions and Coulomb potential are varied at fixed bound- 
ary conditions, while the parameters of the system remain constant. For varying the system 
parameters over the iteration process one needs a proof of existence and uniqueness of solution, 
which is now absent. 

Another approach suggested in [U [6] finds the electrostatic potential f/(r) by means of a 
linear integral equation of the second kind, to which the Poisson equation is artificially reduced. 
In so doing, the necessity of special measures to maintain compatibility of the Poisson boundary 
conditions with charge distribution obtained at each iteration is allegedly eliminated (see [6J, 
Appendix). However it is clear that the exact solution of integral equation, which is equivalent 
to the original differential one, does not exist either, if some of the compatibility conditions 
fail to hold. Apparently for this reason, the said method, even if successful, was complemented 
at each iteration by some unclear renormalization of electron density obtained by means of 
the Schrodinger equation to fulfil the neutrality condition before solving the Poisson equation 
[Ij. One more drawback of this approach appears to be, in fact, the conservation of screening 
term in the region without electrons but with nonzero potential, as it is the case in the surface 
problems. 

In this work we realize the iteration algorithm suggested earlier and described in [:2j, which 
is applicable to solving the Kohn-Sham equations for the equilibrium inhomogeneous electron 
gas and satisfies the boundary conditions for the Poisson equation at any shape of electron 
distribution induced by the effective potential of the previous approximation (section 2). The 
algorithm peculiarities caused by discontinuity of quasiclassical expression for the induced elec- 
tron density with exchange-correlation potential taken into account are pointed out (section 



2.2). The convergence of the algorithm is demonstrated by example of self-consistent calcula- 



tions in a model with homogeneous compensating background while studying such problems 
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as: 1) self- induced bound states near the surface of barrier-bounded electron gas (section 3.1), 
2) parallel-plate capacitor with real distribution of screening charge accounting for Friedel os- 
cillations of electron density (section [3^ ), 3) work function and surface energy of simple metals 
(section |4]). The example of a metal-semiconductor tunneling contact with Schottky barrier is 
used to demonstrate the method's applicability to the systems where the current flow is possible 
(section |5]) . Some of the results were preliminarily presented in and reported to the Russian 
conferences on semiconductor physics in 2001 - 2007. 



2 The base equations and iteration algorithm 
2.1 General formulation 

In the Kohn-Sham approximation, the gas of interacting electrons at temperature T = is 
described by the set of equations for single-particle wave functions 

^V'^eir) + {e- f/efr(r))^,(r) = (2.1) 

with the effective potential 

f/eff(r) = [/(r) + f/.,(r). (2.2) 

Here e is the energy eigenvalue of a single-particle state, and U is Coulomb potential energy of 
the electron described by the Poisson equation 

V^f/ = 47r(A^+(r) -iV(r)), (2.3) 

where iV+(r) is positive charge density, N{r) is electron density. Everywhere, unless otherwise 
specified, the atomic units \e\ = m = h = 1 are used based on charge e and mass m of the 
free electron, when the unit of distance is Bohr radius a^, and the unit of energy is Hartree 
Ha = e'^/ttB- 



The exchange-correlation potential in equation (2.2) is taken in the local approximation of 
density functional 



dN 



(2.4) 

N=N{r) 



where Exc = £x + £c is the sum of exchange and correlation energy per particle. 
The electron density is expressed through wave functions as 

N{r) = 2 j S){£}|^,(r)|% (2.5) 

£<eF 

where factor 2 takes into account the spin degeneracy of single-particle states, and the integra- 
tion on differential spectral measure D {e} of Hamiltonian H = — ^V^ + f/eff(r) is made over 
all occupied states with energy e not exceeding the Fermi energy sp. 

As pointed out in the Introduction, the boundary conditions for the Poisson equation impose 
certain requirements on the spatial distribution of electron density, and the result of integra- 



tion in Eq. (2.5) does not necessarily comply with them until the self-consistency is reached. 



However, if we conceive of the total electron density as the sum 

iV(r) = iViM(f/(r)) + iVq,(r), (2.6) 
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where the induced density A^md depends exphcitly on the unknown potential, the incompatibihty 
problem of boundary conditions with the right-hand side of Poisson equation does not arise 
because the long-range Coulomb interaction is replaced with the screened interaction of finite 
range. In the Kohn-Sham theory, the electron density and effective potential can be coupled 
with an approximate local relation similar to that used for the classical ideal equilibrium Fermi 
gas 

^'^'(^F-f/eff(r))^/^ (2.7) 



ind I 



37r2 



Formula (2.7) can be obtained from the solutions of Schrodinger equation (2.1) for wave func- 



tions of continuous spectrum in quasiclassical approximation by substituting them into Eq. 



(2.5) and averaging the resulting electron density over a scale larger than the Fermi wavelength 
Ap. The function iVqu(r) with A^ind(r) of Eq. (2.7) is a quantum correction to the quasiclassical 



electron distribution. Note that the representation of total density as in Eq. (2.6) is always 



exact regardless of using the approximate formula (2.7) for A^ind(r) 
In expression (2.2), U^c 



Ux + Uc- The exchange potential for electron gas in the local 



approximation is given by the well-known expression 



ox 1/3 



27V J 



2/3 



(2. 



The correlation potential, according to Eq. (2.4), is defined through the correlation energy per 
particle Ec- We assume the latter as 



0.44 



rs + 11.5 

Here rg is the local Wigner-Seitz radius defined as 

47rrf (r)/3 = A^-^(r). 



(2.9) 



(2.10) 



For the following analysis it is convenient to introduce Rs = rs(oo) as the value of this parameter 
at infinity, where we always take the local neutrality condition lim|r|_»oo(^ ~ ^+) — > to be 
fulfilled. 



Having calculated the derivative (2.4) with correlation energy given by Eq. (2.9), one can 
write the correlation potential as 



22 8r,(r) + 69 
'75[2r,(r) + 23]^ 



(2.11) 



The formula (2.9) for correlation energy is a well-known Wigner expression with a different 
parameter 11.5 instead of 7.8 in denominator. The reasons for using the modified Wigner 
formula will be discussed in more detail elsewhere. Here we only note that such a choice of 
Ec, on the one hand, describes better the Rs dependence of cohesive energy of simple metals 
represented in Table 8 of the book [8]. On the other hand, the value Rg = 5.63 accepted in the 
literature for Cs turns out less than the critical one Rsc = 5.64, which characterizes the stability 
of homogeneous electron gas in the jellium model with respect to spatially inhomogeneous 
perturbations. The impossibility to use the formula (2.7) for induced charge in the case of 
Rs^ Rsc would make the suggested algorithm inapplicable, so the choice of correlation energy in 



the form (2.9) provides the existence of cesium metal in the framework of calculation technique 
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under development. There are also other expressions suggested in the literature for correlation 
energy in the local density approximation LDA (see, for example, discussion in [HE]), but for 
the purposes of this work aimed at demonstrating in detail the application of the new algorithm 
to a number of classical problems, the specific form of correlation potential is not so important, 
because the relevant modifications to the expression for Uc lead to only quantitative changes in 
calculation results, and do not create fundamental difficulties for realization of iteration scheme 
suggested. 



The separation of induced charge brings the Poisson equation (2.3) to the form 



V^U + 47riVi,d(f/) = 47r(iV+(r) - iVq,(r)). 



(2.12) 



The quantum correction to electron density Ai"qu(r) introduced here by means of Eqs. (2.5)-(2.7) 
is defined as 

2 J S)H|^,(r)|'-iVi,d(r). (2.13) 

e<eF 



iVqu(r) 



The algorithm of consecutive steps i = 0, 1,2, .. for iterative solving of equations (2.1) - (2.5) 



with equation (2.3) changed to (2.12) can be now formulated as follows: 



1. Neglect the small-scale quantum variations of electron density 

= 0- 

qu ^ ) 

2. Solve the nonlinear Poisson equation with a known right-hand side 

V'U\r) + AnNUU\ iV^J = An (iV+(r) - iV;,(r)) ; 

3. Find the total electron density self-consistent with U, and effective potential 

N^ir) = NUr) + Knir); f/,V(r) = U\r) + U., (iV:(r)) ; 

4. Find the eigenfunctions of the single-particle Kohn-Sham Hamiltonian 

Jv2^:(r) + {e- [/^^(r)) v^^(r) = 0; 



(2.14) 



(2.15) 



(2.16) 



(2.17) 



5. Find anew the total electron density, using the eigenvalues and eigenfunctions obtained 



iV'(r) = 2 J D{£}|vl/:(r) 

£<eF 



(2.18) 



6. Calculate the quantum correction for the next iteration step by means of the new total 
electron density 



(2.19) 



7. Return to solving the Poisson equation with the new right-hand side. 
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Index i here denotes the iteration number. It is of importance that the Coulomb potential 
f/*(r) and the induced electron density N^^^{r) are found self-consistently at a fixed quantum 



density A''* (r) in the process of solving the Poisson equation (2.15). 



The convergence criterion for iteration process is the value of 5, which limits the maximal 
deviation of the electron density N^{r) obtained after solving the Poisson equation and self- 
consistent with the potential U^{r), from the density A^*(r) found by means of the Schrodinger 
equation, 

max \N'^{r) - N\v)\/N{oo) < S. (2.20) 

r ' ' 

It is assumed here that electron density at infinity is determined by the local neutrality condition 
and does not change through the iteration process. 



In the course of solving the self-consistent Poisson equation (2.15), the induced density 



ind 



iVind(f/\iV^ 



should be calculated as an implicit function of unknown potential W at a known quantum 
correction A*"* according to formulae (2.2) and (2.7). It is implicit because the exchange- 



correlation potential in Eq. (2.7) depends on the tota' 



electron density A^j 



ind 



One should note that representing the induced charge as Eq. (2.7) requires the condition 



qu- 



ind 



dU 



< 0. 



(2.2i: 



to be fulfilled. 



Otherwise, the square of the linear screening length in the equation (2.15) linearized for 



U turns negative, and the response to small spatial perturbations of electron density is no 
more localized due to oscillating behavior of solutions of the self-consistent Poisson equation. 



The inequality (2.21) is generally violated in the regions where electron density is small and, 
accordingly, the relative role of exchange-correlation potential increases. In such regions of 
instability, the induced electron density cannot be represented with Eq. (2.7), so we take 



A^*j^^(r) = 0. This situation occurs, for example, near the surface of electron systems at the 
distribution tails where the electron density is close to zero. 

To conclude this section, there are some remarks on terminology used further on. The re- 



sult of self-consistent calculation by means of the complete set of equations (2.14)-(2.19) with 



expression (2.7) for the induced electron density, exchange potential as (2.8 



) and correlation 



potential (2.11) is called the exact or the self-consistent solution. Solving the same set of equa- 
tions, but ignoring f/xc, one obtains the self-consistent solution in the Hartree approximation. 
The self-consistent solution of the Poisson equation (2.12) with the induced density (2.7), ig- 
noring Uxc, and at A^qu = is the Thomas-Fermi approximation. If the exchange potential 
f/x is taken into account in the latter scheme, it becomes the Thomas-Fermi-Dirac approxima- 
tion. Using the total effective potential in the local density approximation with the induced 



charge formula (2.7) is sometimes called the Thomas- Fermi-Dirac-Gombas approximation (see, 
for example, [Qj, Sec. 11). 



2.2 Peculiarities of the algorithm application 

The calculation technique using the induced electron density is clarified more specifically in 
further sections, where the suggested iteration algorithm is realized in calculating the quantum 
correction to the capacitance of barrier structures and in finding the work function and surface 
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energy of simple metals in the model with homogeneous positive background. At the same 
time, the homogeneity of background allows to get still further in the analytical transformation 
of the general equations. Let us consider the semiinfinite electron gas assuming that the ion 
charge forms a homogeneous positive background of density A^_|_ given by 



N. 



0, 



z > z. 

Z < Za 



+ 1 



(2.22) 



In this model, the effective potential in the Schrodinger equation (2.17) depends solely on the 
coordinate z. Such system can represent, for example, a metal surface with self-consistent 
barrier or a semiconductor structure bounded by a potential barrier. 



The expression for induced charge in the form (2.7) is asymptotically correct (with respect 
to quasiclassicality parameter, which is always small in the systems in question at z ^ oo) 
in describing the self-consistent reaction of electrons to the change in the long-wave part of 
Coulomb potential, which is what is necessary to eliminate the infinite range of direct Coulomb 
interaction of electrons. However, as noted before, it fails in the regions of small electron 



density where the condition (2.21) does not hold. In view of the above, the self-consistent 



Poisson equation can be written as 



d^u' 



dC 



+ c„ni,d(«\n;j = c„(^(C - C+) - <u(C)), C > C 

2 =c„(0(C-C+)-ri'(C)), c<c 



(2.23) 



where 



= (S/Stt) {A/9it)^^^ Rs, and 6{() is the Heavyside theta function. We introduced 
the dimensionless variables ( = kpz, = kpz+, n{Q = N{z)/N+, u{Q = U{z)/ep, where 
kp = (37r^A''_|_)^/^, Ep = kp/2. The boundary conditions must force the electric field du/dC, to 



be zero at infinity, to ensure the potential finiteness, m(oo) 
electron gas, n(oo) = 1. 



const, and local neutrality of the 



Eq. (2.23) assumes the existence of a single critical point Q that separates the region with 



the induced density n^ad expressed as (2.7), from the region, where such expression is invalid 



according to criterion (2.21). When calculating the systems that meet the condition < Rsc, 
it is always possible to introduce nind in the region < ^ < oo. The density riq^ turns out 
definite at the same region. Actually the problem is always solved on the finite interval of real 
axis, < ( < Cmax- The value C+ > denotes the possible shift of the positive background 
starting point from zero. Values C+ and Cmax are assumed to be taken so large that the transfer 
of boundary conditions to the points C = and Cmax from ( = —oo and C = oo does not affect 
the result within the chosen accuracy. This condition is easily checked by calculations at several 
consecutively increasing values of C+ and Cmax- 

At the point Cc the right-hand and left-hand solutions are matched together so that potential 
u is continuous along with its first derivative. The critical point Cc is defined as the point where 
inequality(2.21 ) reverses its sign, that is, the derivative dn[^i^/du turns positive. As long as with 



reversing sign the derivative of multivalued implicit function nmd{u) passes through infinity, it 
is more convenient, in the manner of the density functional approach, to analyze the inverse 
single- valued function u{ni^d), whose derivative reverses its sign, passing through zero. This 



function can be found from formula (2.7) and in the dimensionless variables introduced has the 
form 



/i 



n 



2/3 
ind 



[nind + ^qu) 



(2.24) 



7 



where u^c = U^c/^py /i = 1 — Mxc(oo) is the dimensionless Fermi energy of electrons 
accounting for the exchange-correlation interaction. 

Unlike the analysis of validity conditions for the formula of induced density in homogeneous 



electron gas carried out in [lOj, which provides, along with (2.9), the estimate of critical density 
as Rgc = 5.64, in present case the dependence rii^^^u) is also affected by the quantum correction 
to the density. This influence is caused by the exchange-correlation potential in the local 



n. 



qu 



density approximation being dependent on the total electron density. For that reason, the 
analysis should be carried out all over again on the base of the expression 



du 



dn-. 



ind 



-1/3 
ind 



d 



dn 



(^ind + ^qu)] 



(2.25) 



ind 



The result needed for qualitative understanding can be derived analytically, the correlation 
energy being neglected. Its consideration does not change the general picture, for even at the 
least acceptable value of density of the homogeneous electron gas at the infinity, corresponding 
to the highest acceptable value R^c = 5.64, the ratio Uc{Rsc) / Ux{Rsc) = 0.26, and the derivative 
ratio is even less, {duc/ dn) / {du^/ dn) ~ 0.07. With a rise in density, these estimates still improve 
in the sense favorable to our purposes. At the same time, in the course of numerical realization 
of the algorithm the condition du/dni-ad < is easily checked without any simplification by the 



numerically known right-hand side of formula (2.25) 



All the above considered, the inequality that determines the validity range of quasiclassical 



expression (2.7) is conveniently written as 



2 -1/3 
Q^ind 



2 Cx-Rg 



(nind + r2qu)"'/'<0, Cx = 2(2/37r2)'/' 



(2.26) 



The derivative of the exchange potential (2.8) is expressed here in our dimensionless variables. 



On the basis of (2.26) one can obtain the equation for critical values of total density n 



^ind + ^qu, which corrcspoud to changing in sign of the derivative du/dni^d, 



n. 



- (3n, + /3nq, = 0, /3 = (c^i^s) '/S- 



(2.27) 



Its solutions should satisfy the condition of non-negative total density, n > 0, and induced 
density, riind > 0, by virtue of their definition by formulae (2.5) and (2.7). The value of the 
quantum correction be of any sign, because it is defined by formula (12. 61) as a difference 



of two expressions for electron density found from the exact and approximate solutions of the 
Schrodinger equation. 



The roots of equation (2.27) are given by 



^cl,2 



2 



It is clear that with riqu < only one root 



1 ± 



n 



qu 



(2.28) 



2 



1 + 



n 



qu I 1 



(2.29) 



satisfies the stated condition of non-negativity ric and n*-'^^ 



ind 



rip 



riqu. In the case of riqu = 0, 



formula (2.29) gives the value of n^. for the homogeneous electron gas allowing for only the 
exchange potential, which after switching to atomic units corresponds to R^c = 6.02. At 
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^qu < 0, the critical density is even higher. Taking into account the correlation potential also 
increases the estimate for nc, adding a positive term to the left-hand side of inequality (2.26). 
In the case of a homogeneous gas with correlation energy as (2.9), this shifts the critical value 
towards R^c = 5.64. 

While negative values of n^a do not change qualitatively the dependence of induced density 
riind on the average Coulomb potential u compared to the homogeneous case, leaving only one 
value of critical density, there is quite another picture at n^^ > 0. Both roots given by formula 
(2.28) have physical meaning and ensure positivity of both total density ric, which is evident, 
and the induced one Uc — riqu, which can be easily checked. It can also be seen that there is the 
limitation n^^ < (3/4 on the value of quantum correction to provide the existence of the roots. 
If rzqu exceeds this limit, the Coulomb potential u becomes a steadily decreasing function of 
nind in the total range of physically acceptable values, and the applicability condition (2.21 ) for 
the induced density formula (2.7) holds at any values of u and n. 

It should be emphasized that the conditions for ric caused by inequality (2.26) do not 
coincide with the condition of non-negativity of the radicand in formula (2.7), coming into 
action earlier. It is due to this circumstance caused by negativity of u^c and its dependence on 
n that formula (2.7) breaks down earlier with growing potential u, than rii^d turns to zero. So 
at the increasing u or decreasing n, the induced charge drops to zero by a jump from some finite 
value in the course of iterations, creating discontinuities of the first kind in the distribution 
of total electron density and the effective potential, as can be seen in Fig. [TJ After the self- 
consistency is achieved, this discontinuity naturally turns out to be less than accuracy specified, 
i. e. the loop termination criterion (2.20). Nevertheless, the treatment of this discontinuity 
when solving the Schrodinger and Poisson equations turned out an essential factor touching 
upon the effectiveness of the whole algorithm, what will be discussed further. 

The estimate of the value of discontinuities in the density and effective potential can be made 
accurately, neglecting rzqu. It follows from the results of such analysis that the value of density 
discontinuity equals ?t,|^^, taking on the value (i?s/-Rsc)^ at the initial cycle of iterations. One 
can see that with a growing non-ideality of electron gas and Rg closing to Rsc the dimensionless 
density jump tends to unity and becomes a more and more strong perturbation. 



The discontinuity of effective potential equals u^c ny/^ , because the Coulomb potential is 

continuous everywhere. In the atomic units, this discontinuity has an universal value dependent 
only on the expression chosen for correlation energy and ranges from 0.1 Ha with just the 
exchange potential taken into account to O.lSHa with the correlation potential added according 
to the standard Wigner formula [2]. In the dimensionless units, however, the discontinuity 
of effective potential also grows considerably with the increase of electron non-ideality due to 
decreasing e^. 



Now let us write down the algorithm of iterative solution for Eqs. (2.14)-(2.19) as applied to 
the case of semi-infinite many-electron systems with homogeneous positive background consid- 
ered in this section, using the dimensionless variables introduced and duly commenting every 
step taken. 

The initial value of the quantum correction to electron density and the total density itself 
are taken to be zero, and the initial position of critical point is taken beyond the boundary of 
positive background: 



n 



qu 



(C) = 0, n°(C) = 0, C° = C+ - 1. 



(2.30) 



Solving the Poisson equation (2.23), we find the Coulomb potential satisfying the specified 
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boundary conditions and a new total electron density n*, self-consistent with it, 



< = <d + <u, C>C, < = n\C<C, ^ = 0,1,... (2.31) 

In the process, a check for condition dnl^^/du^ < can show that critical point should be 
relocated to a new position Cc^^- In the presence of functions nl{() and m*(C) being numerically 



known, the simplest way of this check is to apply the formula (2.25), where nl^^ is replaced 



with its expression in terms of potential and definition (2.7) is rewritten in dimensionless 
variables: 



ri: 



{fi ~ - u^,{nl)f\ (2.32) 



ind 



If critical point is shifted to the right, one should redefine n*, replacing Q with Q'^ in the 



procedure (2.31). Shift of the critical point to the left leaves unchanged, though the new 
position of critical point defines the region of induced density formation after the Schrodinger 
equation has been solved. 

The case of the initial cycle, i = 0, stands somewhat apart because the direct numerical 
solution of the second-order differential equation can be bypassed. Taking = 0, one should 



begin the solution with Eq. (2.23). In the absence of the quantum correction to density, 
when the induced charge, nonlinear in u, does not depend explicitly on the coordinate (, the 
first integral of this equation can be found in an explicit form. The solution is then found as 
numerical quadrature of the first integral. In that solving one determines the Coulomb potential 
u, the electron distribution self-consistent therewith and reduced here just to n^^^, and the 
critical point position If the system as a whole is sufficiently described by quasiclassical 
approximation, such a solution may have a practical significance (see, for example, formulae 
(35.70)-(35.73) in [H] and the calculation of tunneling current through the Schottky barrier in 

my 

With the solution of the Poisson equation known, one can form the effective potential 
^eff ~ ^* + ^xc(^s) th^t enters into the Schrodinger equation. In virtue of translation invariance 
of the system under study in [x, y) plane, the total wave function can be written as 

^fc,k||(C,r||) = 7^e^'^Nr|ic^,(^), (2.33) 

where vector r|| lies in the (x, y) plane, and constant C is defined by the normalization require- 
ment for function ip^. The Schrodinger equation for ipl at the i-th iteration cycle is conveniently 
written as 

dC 

where the dimensionless quantum number k is normalized to kp, and the scalar potential is 
calibrated by the condition ■u(oo) = 0. An eigenvalue of k specifies the behavior of continuous 
spectrum eigenfunctions normalized to 6{k — k') according to [2] (see formulae (24)-(25)) with 
the asymptotic form at infinity 

i^iiO^ ylMK + lk), C^oo. (2.35) 



+ {k' - <ff(C) + «xc(oo)) ^liC) = 0, (2.34) 



The numerical solution ipkiC) the Cauchy problem for the homogeneous equation (2.34) 



found accurate within an arbitrary factor, is normalized correctly, if its asymptote at C ^ oo 
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is as Eq. (2.35). This is achieved multiplying the obtained solution by a constant A defined by 

^'(v;^(Coo)+^^'(Coo)) =2/7r. (2.36) 

Here (oo < (max is a sufficiently remote point where potential Meg (Coo) can be considered as 
become a constant. The phase jk can be found from the relation 



7fc = arctan 



mo 



(fcCoo - ^vr) 



(2.37) 



where / = [k(c>o/7f] is an integral part of the ratio. With such a definition, the phase is invariant 
with respect to the choice of Coo point. 

To conclude the i-th iteration cycle, the total electron density is calculated. 



n^(C) = 3 f dk{l - k') lijliOl 
Jo 



as well as the quantum correction to the density for the next cycle 



(2.38) 



(2.39) 



with the definition (2.32) taken into account. 



The Cauchy problem for the Schrodinger equation was solved numerically using the fourth- 
order implicit technique that is often called the Numerov's method. The Numerov's algorithm 
does not include directly the value of first derivative of the sought solution specified at the initial 
point. Instead of this, one should know the value of the sought function at the first point from 
the boundary. To calculate this value with accuracy matching that of the Numerov's method, 
we represented the function as a second-order expansion in Taylor series. The necessary value 
of the second derivative at the boundary was found by means of the Schrodinger equation 
itself. A similar technique was applied to pass the discontinuities of effective potential. It has 
been checked within the chosen accuracy that the result of such a procedure is in accord with 
Runge-Kutta method commonly used to move one step away from the point where the Cauchy 
boundary conditions are specified. 

The Poisson equation was solved by the relaxation method fTT], its idea being to solve the 
matching nonstationary equation 



du 
di 



+ Cn^ind(M,'^qu) = -Cn[0{C - C+) 



n, 



qu 



(C)] 



(2.40) 



instead of the stationary problem, until the solution stops changing in time t within the accuracy 
specified. 



From the formal point of view, Eq.(2.40) is a nonlinear diffusion equation, a lot of numerical 



schemes having been developed for its solution. In present work we applied a scheme, in which. 



after discretization (2.40) in time with step r and in space with step h, one specifies an initial 



potential distribution satisfying the boundary conditions, and time evolution of the solution 
is monitored by solving the boundary value problem at each time point. As this takes place, 
one need not to obtain with a good accuracy the spatial distribution of potential u at each 
"time step", searching only for an accurate stationary solution at long times when change in u 
from step to step becomes small. The time discretization was made with an implicit first-order 
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scheme, so, when forming the boundary value problem at a new time point, the nonlinear term 
in the equation was estimated via the induced density ^ind(M, ''^qu) linearized in the increment 
of potential at the time step r. In the course of solving Eq. (2.40), the choice of r value was 
dictated, on the one hand, by considerations of numerical stability, and on the other hand, 
by those of computation speed (the time of reaching the stationary solution). In line with 
recommendations given in the literature, the initial time step was chosen according to relation 
r//i2 = 0.1. 

Since the linearization of dependence nind(^) was used at the movement from one time point 
to the next, the necessity could arise in the course of solving the Poisson equation to redefine the 
critical point to avoid the appearance of a region where dniad{u,nq)/du > 0. This was carried 
out similarly to the procedure described after formula (2.32). As noted above, the magnitude 
of ordinary jump of the electron density increases with Rg and leads to more and more violent 
perturbations of numerical solution. It was found (see details in Sec. Ill) that at Rg < Rsc the 



convergence of solution of the equation set (2.14) - (2.19) deteriorates. It appears though to be 
rather a consequence of some specific numerical algorithm realization, than a drawback of the 
iteration scheme suggested. 



3 Properties of electron gas in the barrier structures 

In this section, we apply our method of solving the Kohn-Sham equations to finding the self- 
consistent electron density and effective potential in a barrier structure that consists of two 
conductors with degenerate electron gas separated by an insulating layer with dielectric per- 
meability Hd- We assume that z-axis is perpendicular to the layer plane, the conductors are 
positioned at z < —d/2 and z > d/2 respectively, and the insulating layer lies at \z\ < d/2. It 
is assumed further on that the neutralizing background in the conductors is homogeneous, and 
the barrier is not impurity-doped, that is. 





\z\ 


[N+, 


\z\ 



The height of potential barrier formed by the insulator is taken to be infinite, so that wave 
functions turn to zero at the conductor-insulator interface. 



3.1 Self-induced energy levels in semiconductors 

In the case of infinitely high potential barrier, the electron states localized in z and the size- 
quantized two-dimensional subbands can emerge even if inhomogeneities of positive background 
N+{z) and applied voltage are absent. The existence of such states was suggested for the first 
time in [12], where the discrete energy levels in question were called "the plasma levels". The 
appearance of localized states was attributed to the electron wave functions turning to zero 
at the barrier interface. Due to continuity of wave functions, the electron density close to the 
barrier is small and inadequate to provide the local neutrality. The uncompensated positive 
charge near the interface creates a self-consistent potential well, wherein bound states can exist. 
According to calculation carried out in [12] for the case of degenerate electron gas, the discrete 
levels should exist at Rg < 2.87, while at Rg < 0.15 the second bound state appears in the 
near-barrier well. 
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Independently of [I2] , the idea of depletion region near the barrier was also suggested in [13] , 
where the energy of localized states and the depth of self-induced potential well were calculated 
as functions of Rg. However, in these works either the self-consistency has not been reached, 
as in [121 , or some model potential was used in the Schrodinger equation instead of the exact 
solution of Poisson equation ^13j. The exchange-correlation interaction of free charge carriers 
was neglected in both papers. The self-induced levels are quite shallow with energies less than 
0.1 of the potential well depth (~ 2 meV for a semiconductor with n-GaAs parameters), thus 
being very close to the states of continuous spectrum. A small change in the potential well 
shape can make them disappear. In this context, it was interesting to find out if these states 
remain in the self-consistent solution. 



The Kohn-Sham set of equations was solved according to iteration scheme (2.14) - (2.19). 
The total dimensionless electron density at z > d/2 obtained after solving the Schrodinger 
equation on the i-th iteration was calculated by formula 

n\0 = 3 / dk{l - k'MiOl' + ^ E I^KOr(l - e]), (3.2) 
Jo ^ ^ 

where ( = kp{z — d/2), j is the number of a discrete level, ej are discrete energy levels, ipk and 
ijjj are the respective wave functions of continuous and discrete parts of the spectrum. The first 



term in(3.2) corresponds to the density of continuous-spectrum electrons, and the second one 
is the electron density in the localized states. 

Since the considered problem with an infinitely high barrier is close to that of potential 
scattering of s-wave by a short-range potential, one could search for energy levels in the potential 
well by means of the Levinson's theorem (see, for example, [2]) that expresses the number m of 



levels in a well in terms of the asymptotic phase (2.35) of continuous-spectrum wave functions 
at infinity by relation limfc^o 7fc = '^tt. Although this formula is used, for example, in [15] for a 
theoretical analysis of barrier structure, it is not convenient in numerical calculations because 
of redundant computation of wave functions for empty states. Instead of this, it suffices to find 



a solution of the Schrodinger equation (2.34) at A; = and calculate the number of zeros within 
the interval < C < 00. 

However, when finding the energy and the wave function of a bound state we applied the 
trigonometric sweep method, which eliminates the well-known numerical instability of solu- 
tion decreasing in the classically forbidden region, caused by the second solution that grows 
exponentially [^. In this method, the wave function is written as 

^(C) = a(C)sinr?(C), ^ = a(C) cosr/(C), (3.3) 

so one can easily show that the total phase r]{() introduced here has the increment vr at zeros 
of the function ip{C) where the logarithmic derivative of the latter turns to infinity. There are 
two first-order equations for the phase and the amplitude 

dTj 

— = {e- Meff(C) + Meflf(oo)) siu^ T] + COS^ T] (3.4) 

^ = sin 2r] (1 + Ucs{C) - Ucs{oo) - e) , (3.5) 



which are derived from (2.34) taking (3.3) into account. The equation for ri{C,) is solved with 



initial condition 7^(0) = 0, while the initial condition for the amplitude can be taken at random 
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and specified definitely by the wave function normalization to 1. The level energy Sj < 
corresponds to the limit value of phase ^{oo) = jvr, (j = 1, 2, ..). So to calculate the number of 



levels in a well, one should only solve the equation (3.4) for two energy values corresponding to 
the potential well bottom and the edge of continuous spectrum e = 0, and find the tt multiplicity 
of the difference between the two phases at infinity. 




Figure 1: Electron density obtained at various cycles of iteration process. Number of curve 
corresponds to the number of iteration after which it was calculated. In this section, coordinate 
z in the figures is counted from the point d/2. 



Boundary conditions for the Poisson equation were taken as 



du 



C=o 



du 
dC 



(3.6) 



Note that such boundary conditions along with Fermi level specified by the requirement of local 
neutrality aX ( oo and the choice of u{oo) = define a unique numerical solution of the 



Poisson equation (2.23) and ensure its tending to zero at infinity. 



FigjT] demonstrates the electron density distributions obtained by various number of iter- 
ations. One can easily see that, at the inintial iterations, the spatial dependence of electron 
density has discontinuities at the boundary of validity of quasiclassical expression for induced 



charge in the form (2.7). The discontinuities disappear in the course of iteration process, and 
final self-consistent density is a smooth function of coordinate. It is also essential to note that 
iteration convergence criteria specified were usually reached at 10-14 cycles (depending on Rg). 
The continued calculations did not provide a noticeable increase in accuracy at chosen dis- 
cretization parameters, but neither disrupted the solution, as can be seen from Fig. 1, which 
testifies to the stability of the algorithm applied. 

Fig. [2] demonstrates the coordinate dependence of the self-consistent potential and electron 
density. As can be seen from Fig. [2]A., taking into account the t/xc modifies shape and depth 
of the potential well and, correspondingly, the energy of localized state (see curves 1 and 2 in 
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Figure 2: Self-consistent potential and electron density at Rg = 0.4. Fig. 2A demonstrates the 
self-consistent potential with account of the exchange-correlation interaction (solid line) and 
that in the Hartree approximation (dashed line). Fig. 2B shows the total self-consistent electron 
density (curve 1) and density of electrons localized in the size- quantized subband (curve 2). 

Fig. |3]). The total self-consistent electron density in Fig. |2]B (curve 1) is the sum of density 
of electrons in the size-quantized subband (curve 2) and that of electrons in the continuous 
spectrum states. Localized electrons are concentrated mainly close to the potential barrier, 
though the localization pattern depends on Rs. 

The Rs dependencies of self-consistent well depth and energy of the state localized therein 
are presented in Figs. [SjA. and |3]B respectively. Curves 1 and 2 in both figures are results 
obtained in present work, while curves 3 and 4 are results presented in papers [12] and [13]. 
One can see that calculated level energy and its Rs dependencies presented in [12] (curve 3B) 
differ both quantitatively and qualitatively from our self-consistent calculations. Furthermore, 
we did not found two localized levels in the potential well in the absence of external electric 
field at any values of Rs > 0.05. The second level with a dimensionless cohesive energy about 
0.0016 was found only at Rg = 0.005 that is not, as a rule, realized physically in real structures. 
Nevertheless, one should note that the idea of a self-consistent potential well with bound states 
close to infinitely high barrier, put forward in [i2j, proved to be true. 

The results reported in [13] for the potential well depth are in rather good agreement with 
the self-consistent calculation of the present work in the Hartree approximation (see. curves 2 
and 4 in Fig. |3]A). However, the energy of bound state (curve 4 in Fig. [3|3), calculated in [13] 
differs from the accurate solution both in magnitude and the Rs dependence, particularly at 
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Figure 3: Potential well depth (A) and bound-state energy (B) as functions of Rg. Curves 
1 were calculated with account of the exchange-correlation interaction, curves 2 are results 
of self-consistent calculation in the Hartree approximation. Curves 3 and 4 represent results 
presented in |T2] and [13] respectively. 

small values of this parameter. At the same time, the difference cannot be attributable to the 
exchange-correlation interaction neglected in [13], as it is evident from the relative closeness of 
our curves 1 and 2 in Fig. [3|3 to each other. 

It may be safely supposed that considerable error in the self-consistent energy levels found 
by method used in [TBj is caused by inaccurate description of the potential well shape by a 
three-parameter model potential. 

Taking into account the finite barrier height leads to a change in boundary condition for 
electron wave functions at the insulator/conductor interface, and that, as calculation shows, 
leads in its turn to a lowering upper limit of Rg region where localized states are observed. At 
a considerable decrease of barrier height the discrete level disappears entirely. 

Far from potential barrier, kpz ^ 1, the effective potential and electron density manifest the 
Friedel oscillations (see. Fig. |4]). One can see that the oscillation amplitude of effective potential 
with f/xc component is considerably less than that of Coulomb potential. Such suppression 
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Figure 4: Friedel oscillations of effective potential and electron density in the electron gas 
limited by an infinitely high barrier. Curve 1 is effective potential, curve 2 is the Coulomb 
component of effective potential, curve 3 is electron density 



of oscillations can be understood as follows. According to Poisson equation, the oscillating 
component of electron density is maximal at the same points as the Coulomb potential. The 
exchange-correlation contribution into Ues always has a sign opposite to that of the direct 
interaction. So in virtue of the repulsive nature of electron-electron interaction, the potential f/xc 
is negative, and its magnitude grows steadily with electron density according to formulae (2.4), 
(2.8) and (2.11). Consequently, at the points of maximum Coulomb potential the exchange- 
correlation contribution is also maximal in magnitude and negative. 

The similar considerations show that at the minima of the effective potential the Coulomb 
contribution and the exchange-correlation one also summarize in phase opposition. As long 
as the relative role of f/xc grows with decreasing electron density, this leads to the oscillation 
magnitude of effective potential being suppressed more than that of the Coulomb one with 
growing Rg. Note that the influence of t/xc upon the Rg dependence of oscillation magnitude of 
the dimensionless electron density is quite opposite. The reason for this is discussed in more 
detail in Section ID 

In conclusion, one should note that the system considered with filled states of discrete 
and continuous spectrum has much in common with quasi two-dimensional electron gas of 
accumulation layers close to semiconductor-insulator boundary with an external electric field 
applied or even without it. Though one believes [T7j that the problem of self-consistent spectrum 
calculation for size-quantized subbands in the enriched layer is solved in publications 
concerned with methods for calculating such structures still appear [THj. However, as before, 
no proved technique is suggested allowing to avoid the mismatch between the total charge or 
the dipole moment of quasi two-dimensional gas obtained at intermediate steps of iteration and 
boundary conditions to Poisson equation imposed on the potential. To make a solution possible, 
the authors renormalize the electron density found by means of the Schrodinger equation, and 
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also vary the Fermi level position in the bulk of semiconductor. The results obtained in the 
present section of our work make clear that the algorithm suggested allows one to obtain the 
solution in a regular manner within the framework of mathematically complete and accurate 
iteration procedure. 

3.2 The capacitance limitation of a barrier structure 

Up to now it was assumed that external voltage is not applied to barrier structure. In this 
section, we find the self-consistent distribution of electron density and effective potential in the 
case of a dc bias V, and calculate the capacitance of the barrier structure. Apart from the check 
for the algorithm performance, this problem is of independent interest in connection with the 
capacitance limitation of thin-film capacitors in general, posed in [19], as well as with particular 
case of nanocapacitor with a thin-film ferroelectric insulator in view of its application prospects 
in electronics pO]. 

It was shown by Mead [12] that the measured inverse capacitance of the metal-insulator- 
metal junction as a function of insulation layer thickness d does not tend to zero at d 
in contrast to the common formula for geometrical capacitance of a plate capacitor. This 
was explained qualitatively by the finite thickness of the spatial charge layer that screens the 
electric field penetrating the real metal. The quantitative analysis of this situation in the 
framework of the Thomas-Fermi approximation for electrons in the metal was carried out in 
f2T\ and complemented with an account of dielectric permeability of the metal electrode in 
[22] . In the latter work, it was also noted that using the insulator with a high permeability 
Kd makes the effect of non-ideality of electrons more pronounced even at insulator thickness 
much greater than the Thomas-Fermi screening length in the electrodes, and the maximum 
achievable junction capacitance with given electrodes was introduced as the limit at — > oo. 

The formulae for capacitance obtained in the above papers provided a qualitative under- 
standing of its experimental dependence on the material of electrodes and insulator, and also 
allowed one to get a correct order-of-magnitude estimate for the contribution of spatial charge 
region to the capacitance. However it was noted that the effect of the exchange-correlation 
interaction of electrons in metal and that of the charge in the Friedel oscillations of the density 
are still unclear and cannot be accounted for in terms of the Thomas-Fermi approximation. 
At the same time, as early as in work [23] by Bardeen it was actually shown that spatial elec- 
tron redistribution related to the density oscillations caused by the barrier does not keep the 
total neutrality of the system because of the occurring deficit of electron charge. In the case 
of an infinitely high barrier without an external field applied this charge deficit equals (in the 
approximation of non-interacting electrons) 



which coincides with the excess of the positive charge given by formula (26) in [23] (with 

printing error corrected, see. Eq. (2) in [24j). Besides, it follows from the results of previous 
section that self-induced surface level also can capture a considerable electron charge (see. Fig. 
[2^). With a noticeable field dependence of charge in the oscillations and on the near-barrier 
level their contribution to the capacitance may be comparable with that of the charge in the 
Thomas-Fermi screening region. 



oo 




(3.7) 
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All above factors can be estimated by solving self-consistently the Kohn-Sham equations for 
a barrier structure in the jellium model. Let us consider a system of two semi-infinite metals 
separated by the infinite potential barrier of width d and dielectric permeability kj. Such a 
barrier prevents the charge transfer between two parts of the system when voltage is applied. 
The differential electrostatic capacitance per unit area of a wholly neutral system is given in 
terms of the charge Q per unit area of the right-hand side of the structure by formula 

dO d 1"^ 

CiV) = -^ = ^J^ dz{N{z, V) - N^{z)), (3.8) 

where V = ip(—oo) — (p{oo) is voltage across the structure equal to the work done by the electric 
field on a unit positive charge to move it from negative to positive infinity. Here (p{z) is the 
dimensional Coulomb potential. The point Zq given by dD/dz\zQ = is the sign changing point 
of charge density in the region where electric induction D is maximal. As long as electric field 



is zero at both infinities because of the finiteness of potential (fiz), the definition (3.8) provides 
that charges in two parts of the structure are equal in absolute magnitude. 

We assume that the metal electrodes are identical and the positive background distribution 



is described by formula (3.1). To simplify calculation, we also neglect a possible nonzero value 
of the contact potential difference along with related built-in fields and charges in surface states, 
as it was assumed in the model considered in [21] and [22]. These conditions mean the absence 
of charge within the barrier, which allows one to take any point of the insulator as zq including 
the boundary ones zq = ±d/2. If the induction D inside the barrier is specified, the right-hand 
electrode charge equals 

QiV) = -D/Ati. (3.9) 

We represent the voltage across the total structure as a sum of voltages across its three con- 
stituent parts, 

V = [^(-oo) - ipi-d/2)] + E^d + [(^(d/2) - <^(oo)], (3.10) 

where field inside the barrier = D/k^- Let (f±{z,E) be the Poisson equation solutions at 
d/2 < z < oo and ~d/2 > z > — oo with boundary conditions d(p±/dz\z=±d/2 = E, (f±{±oo) = 
respectively. Then, in view of induction continuity at the insulator-metal boundary, the 



equality (3.10), can be rewritten as 

V = -ip4-d/2,D/Kj + Edd + ip+{d/2,D/Kj. (3.11) 

Let Eja = D / be the field in metal at the insulator interface, where the not-equal-to-one 
dielectric polarization takes into account a possible polarizability of internal electron shells 
of ion cores. Since in our case (p_{—d/2,E^) = {p+{d/2, —E^), we get a resulting formula, 
which reduces the calculation of barrier-structure capacitance to the problem for semi-infinite 
system considered in the previous section, 

V = v+{d/2, E^) - -E^) + E^d. (3.12) 



We calculate the capacitance at small voltage when differentiation in Eq. (3.8) can be replaced 



with the ratio of small increments. Due to the additive structure of V in (3.12), it is more 



convenient to deal with the inverse capacitance. Taking (3.9) into account, we get 



C-^^y ^^fMXM^pmzM + 4, A , 2cr' + Ci\ (3.13) 
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where the standard notation Ci, Cd is used at the end for separating contributions of voltage 
across the interfaces (inside the metal electrodes) and the insulator to the total capacitance. 



The explicit definition for Ci and can be easily seen from the structure of expression (3.11), 
if one takes into account that ip+{d/2,0) ^ (see. Fig. [2jA.). It is evident that the maximal 
possible capacity Cmax of the considered structure is achieved in the limit of small electric 
thickness of insulator, limrf/Kd — >• 0, and can be written as 



K 



E 



(3.14) 



We assume further on that = 1, for, as may be required, the numerical value of capacitance 
being found can be easily recalculated using Eq. (3.14) and taking into account the relevant 
redefinition of atomic units. 



The structure of formula (3.14) shows that, for calculating the capacitance limitation of 



barrier structure one should solve twice the problem for semi-infinite metal bounded by an 
infinitely high barrier with boundary conditions as 



du 

dC 



±E, 



C=o 



du 

dC 



0. 



(3.15) 



Here we use the dimensionless variables for electron potential energy and spatial coordinate 
introduced in Section 12. 2[ The dimensionless electric field E is normalized to intrinsic field 
Ec = kpEp/e. Apart from the change in the boundary condition at zero compared to (3.6), no 
modifications were needed in the calculation program. 

The found values of potential at the insulator boundary at two values of electric field were 



used in formula (3.14) for capacitance calculation. For the purpose of accuracy and to verify 



that the chosen field magnitudes are sufficiently small and do not exceed the limits of linear 
voltage dependence of the charge, the potential was calculated for three pairs of values ±E, 
and the slope of the linear dependence obtained was used to calculate the capacitance. 

Fig. [5] demonstrates the Rs dependence of capacitance limitation (denoted just as C) of 
the structure obtained in different approximations. Curves 1 and 2 correspond to the self- 
consistent calculation of capacitance with account of the exchange-correlation interaction and 
in the Hartree approximation. Curve 3 was obtained in the Thomas-Fermi approximation and 
can be represented by a simple explicit formula 



C 



h 



TF 



1 



Stt 8 



1/3 



R. 



1/2 • 



(3.16) 



where /ctf is the inverse Thomas-Fermi length. If, in the Thomas-Fermi approximation, we 



also take into account the exchange-correlation potential in accordance with formula (2.7), we 



get the capacitance formula like (3.16), but with /ctf replaced with renormalized value /ctf 



fcTF/(l + ^durcc/dny/'^ . The corresponding curve is not shown in Fig. [sjbecause the difference 
between these two curves in the considered range of Rg is almost just as small as between 
curves 1 and 2. With an exact calculation, the account of exchange-correlation potential leads 
to the decrease of capacitance, which is probably related to withdrawal of the minimum of 
self-consistent potential well and charge accumulated therein from the barrier (see Fig. 2a). 

The comparison of numerical values of interface capacity Ci = 2C obtained from curve 
1 with its experimental estimations (see, for example, [25j) shows an agreement in the order 
of magnitude. Thus, not taking into account the dielectric permeability of electrode metal, 
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Figure 5: Capacitance of barrier structure at various values of Rs- Curve 1 is self-consistent 
calculation with account of exchange-correlation interaction, curves 2 and 3 are the same in the 
Hartree and Thomas-Fermi approximations respectively. 



which differs noticeably from unity for the noble and transition [27j metals, the range of 
Ci change, according to curve 1 in Fig. [sj is from 275 to 85 fF//i^ at 0.5 < Rs < 3.5. In 
terms of the effective electric thickness dcs = (47rCj)^^ we get the range from 0.035 to 0.11 nm. 
The analysis of current theoretical and experimental data made in [28] led to the conclusion 
in favor of imperfect screening in metal contacts as the main reason for the so-called "dead 
layer" observed experimentally. The latter determines the maximum achievable capacitance in 
structures with ferroelectric insulator for a specified choice of metal for electrodes. 

The general conclusion from the data presented in Fig. |5] is that the Thomas- Fermi approx- 
imation with its exponential distribution of spatial charge does provide the order-of-magnitude 
agreement with complete self-consistent calculation as concerns the capacitance limitation. Nei- 
ther the exchange-correlation interaction of electrons, nor density oscillations or an account of 
the self-induced potential well with a discrete energy level therein prove essential. The level 
that appeared at Rs < 1.6 did not affect the smoothness of the curve C{Rs), in accordance 
with result reported in [29]. The similar absence of drastic changes in the experimental curve 
C{V) at changing number of size-quantized levels in the accumulation layer of n — InAs was 
reported in [HUj . 

In work [15], the preservation of continuous dependence of self-consistent well potential 
on the external electric field regardless of the bound state appearance or disappearance was 
explained by forming of an "orthogonal hole" in the spatial distribution of electrons filling 
the states of continuous spectrum . Our numerical results demonstrate that this phenomenon 
is determined by appearance or disappearance of one more zero in the continuous-spectrum 
wave functions with the appearance or disappearance of a bound state in the potential well. 
As a consequence, filling of a bound state is accompanied, in a sense, by redistribution of 



total electron density among its two components in formula (3.2), rather than contribution 
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from the level is added to the existing one from the continuous spectrum. On the contrary, 
failing to account for (or neglecting) the contribution from a bound state, even quite shallow, 
in the iteration process leads to formation of a shallow and wide potential well with more levels 
appearing at the next iteration cycle, so the iteration process fails to converge. 

4 Work function and surface energy density of metal in 
the model with homogeneous background 

In the model with homogeneous positive background (jellium model), the metal work function 
and surface energy density were calculated in a good many papers. However, the existing 
data reveal sometimes not only quantitative, but also qualitative difference. We believe that 
the reason for such divergence is related to the complete self-consistence being not reached in 
the works under discussion when calculating the distribution of electron density and effective 
potential. Among the detailed and comprehensive works with calculation technique described 
one could mention publications |3], [31] and [32], our results we will compare with. Note that 
the purpose of this part of our work is to obtain an actually self-consistent solution for one of 
the classical problems of the theory of inhomogeneous electron gas and assess the importance 
of self-consistency by comparing the results with those obtained by other methods for the same 
problem statement. We do not aim at improving the jellium model or the density functional 
theory for calculating the surface properties of metals. 

The problem statement here is generally similar to calculation of barrier structure properties 
considered in the previous section. However, in the case of metal the surface barrier is formed 
entirely by the self-consistent field and does not contain any specified seed potential. Besides, in 
metals in ordinary conditions, the values of Rg calculated for charge carriers with free electron 
mass in the medium with dielectric permeability = 1 exceed considerably Rg values typical 
for semiconductor structures and lie in the range 2 < i?s < 5.7. All the above raised the interest 
for a calculation intended to explore the behavior of convergence to the self-consistent solution 
and assess the algorithm stability at Rs values exceeding those considered in Section |3} 



4.1 Work function of simple metals 

Let the positive background occupy the half-space z > z+ where z+ is coordinate value at the 
metal-vacuum boundary. In the calculation, the boundary conditions for the Poisson equation 



(2.23) were taken as du/d(\i^=±oo = 0. Wave functions ipkiC) of single-particle states with 



energies not exceeding the Fermi level were found as Cauchy problem solutions exponentially 



damped at C ^ —oo and normalized according to Eqs. (2.33) and (2.35). The dimensionless 
work function w = W/ £p (see Fig. [6]) equals the potential barrier height counted from the value 
of bulk chemical potential, 

w = Mcfr(-oo) - fi. 

The self-consistent electron density and effective potential obtained after convergence of itera- 
tion process in accordance with calculation scheme described in Section [2] are shown in Fig. [6] 
The curves n and Ues taken from the table 1 in the work [3j are presented by dots in the same 
figure for comparison. One should note that in work [3] the correlation energy was taken as 
standard Wigner formula 
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Figure 6: Comparison of self-consistent calculation for effective potential and electron density 
(solid curves) in present work with similar data from work [3| (points) at Rs = 4. 



rather than expression (2.9) used in present work. For the correct comparison, the result of our 



self-consistent solution presented in Fig. [6jwas calculated using correlation energy (4.1). On 
can see that calculation results in the figure agree rather well, though in work [3] a 12-parameter 
formula for electron distribution was used instead of searching for self-consistency, and those 
parameters were found by some variational procedure (see. Appendix B in [3]). 

However, the convergence of our iteration solution and its good agreement with the generally 
accepted calculation results from the work [3J by Lang and Kohn at a particular value of -Rs are 
not in themselves sufficient for a full assessment of the self-consistency degree achieved. In the 
work [33] within the model of homogeneous background, an universal relation for electrostatic 
potential was obtained, 

U{z^) - f/(oo) = N— ^ Abv, (4.2) 

where e is total energy per one electron in homogeneous electron gas. For further convenience, 
we introduce the dimensionless parameter of self-consistency Abv = ^bnI^y ^^"^ write it as 

2 / 4 \ 

+ 2 - Rl[U^^-e^,], (4.3) 



where the expression in square brackets is left in atomic units. The relationship (4.3) was 
used to verify if the calculated potential is actually self-consistent with the electron density. 
Good agreement between A = m(C+) — u{oo) and Abv proves that a genuine self-consistent 
solution of the Kohn-Sham equations was found (see. P- 1430). The Table [l| presents 

the values of A and parameter Abv corresponding to the results of self-consistent calculation 
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Table 1: The Rs dependence of parameters A and Abv as well as work function W and surface 
energy density S. The number of iterations needed to reach self-consistency is also presented. 
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for specific Rg and values in accordance with formula (2.9). One can see that A and Abv 
match for all Rs < 3.99 to an accuracy better than 0.001. Contrary to work [32] where the 
criterion (4.2) was used explicitly in the process of solution building, this relationship is never 
taken into account in our algorithm. So its satisfaction proves that n and u^s obtained are 
genuine self-consistent solutions of Kohn-Sham equations for assumed approximation for the 
exchange-correlation potential of electrons. 

One should note the sensitivity of the self-consistency criterion to the values of potential at 
the positive background frontier. The degree of agreement between A and Abv presented in 
Table 1 is achieved only if the potential m(C+) value taken is interpolated to the discretization 
half-step 6( towards vacuum region. In the calculation, 6( was chosen within the interval 
0.05^0.005. 

A similar self-consistency test for electron density and potential found by Lang and Kohn 
[3] was carried out in the work [33]. The data presented there in Table 1 show that, even in the 
case of improved (unpublished) results of Lang, the self-consistency criterion (4.2) is fulfilled 
to an accuracy better than the third decimal place only for Rs < 4. With growing Rs the 
agreement is steadily worsening and aX Rs = Q the deviation reaches 0.01. This means that at 
large values of Rg close to the existence limit of real metals, the parametrization scheme chosen 
in work [3] and the method of calculating the parameters are not quite an adequate replacement 
for an accurate solution of the Kohn-Sham set of equations. This fact can be considered as 
an indication of serious difference in the behavior of electron density and potential at the Rg 
values close to critical Rgc, from the jellium model solutions at higher densities of electron gas, 
which was discussed in some measure in Section |2l 

It follows from the Table 1 that in our calculations at /2s — 5 the self-consistency criterion 
is also fulfilled worse than at small values. In the region Rs < 3.99, the iteration process 
was programmed to terminate if the discrepancy defined by formula (2.20) becames less than 
specified value 6 = 10~^. However, at Rs = 4.96 the algorithm of critical point (c autodetection 
described in Section 2^ failed to provide convergence of the Poisson equation solution due 
to the steady movement of critical point Cc from the metal surface into the bulk. To obtain 
a definite result in this case, one had to turn off, after a number of initial iteration cycles. 
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the program module that redefined (c- Nevertheless, after such compulsory choice of critical 
point, the electron density jump in its vicinity did not drop lower than 3 ■ 10~^ in the course 
of iterations. At i?s > 5 the convergence of iteration solution was not achieved at all. It 
appears that the crucial factor here, as concerns the algorithm operation, is that at Rg — >■ -Rsc 
the effective potential (2.2) turns out almost equal to the local exchange-correlation one, and 
the contribution of the long-range Coulomb potential to forming the surface barrier becomes 
actually inessential. 

The Table 1 demonstrates also the calculated Rg dependence of the work function W. In 
Fig. [Tjthis dependence (curve 3) is compared to similar ones presented in works pi] and [32] • 
The values of work function calculated in present work are in rather sensible agreement with 
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Figure 7: Work function in the homogeneous background model. Curves 1 and 2 represent 
results of calculations in [31] and [32]. Curve 3 was calculated in the present work. Circles 4 
show experimental values of work function for polycrystalline specimens of Al, Li, Na, K and 
Rb from Table 3 of review [31]. The smaller of two W values at Rg = 2.07 (Al) was measured 
in polycrystalline Al in the work 



experimental data for alkaline metals (Li {Rg = 3.28), Na {Rg = 3.99), K {Rg = 4.96)) and 
lie closer to measured values than those of curve 1 obtained in [31]. Most likely, the improved 
quantitative agreement between our calculations and experiment is a consequence of using the 



modified Wigner formula (2.9) for correlation energy, because at the same choice of expression 
for Ec the results for the calculated potential are close to each other at least for Rg = 4, as can 
be seen in Fig. |6} 

Our results for the work function at i?^ > 1 are in semi-quantitative agreement with those 
obtained in [32] (Fig. |7( curve 2), and this agreement improves with growing Rg. One can 
see from the figure that both curves have maxima at i?^ ~ 2, however, at small Rg there is a 
qualitative difference. According to [32J, work function has a local minimum at i?s ~ 1, which 
is not confirmed by our calculation. This local minimum in work function is possibly related to 
the absence of self-consistency in the solution of Kohn-Sham equations that has been obtained 
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in work [32] for small values of Rs. 

The method of solution used in [32] and described in detail in work [36] replaces the Poisson 
equation with an integral equation for potential in some finite vicinity of metal surface and 
apphes boundary conditions that do not fix the total charge (see. [3S], Appendix A). Along with 
the problems noted in Section [T] in relation to the replacement of Poisson equation with integral 
one, such boundary conditions cannot provide a self-sustaining consistency of the final result of 
iteration procedure, as it was supposed in [32] (see the text after formula (7)). As mentioned 
in [32], the process failed to converge at all at Rg < 0.5. The reason for this lies probably in 
the increasing role of Coulomb potential compared to that of the exchange-correlation one and, 
therefore, the increasing sensitivity to self-consistency with a rise in electron density. 

We note once more that the work [32] and ours are identical in the physical problem state- 
ment and in the formulation of initial equations except for different correlation energy formula 
used. At the same time, we faced no crucial difficulties with convergence at the region of small 
Rg and were limited only by growing requirements for computer memory and increased compu- 
tation time caused by widening vacuum region with a non-zero electron density. As concerns 
the calculation problems in the region of large Rg, they require a further analysis and possibly 
relate in part to the use of local density approximation for the exchange-correlation potential. 
On the one hand, a ground for such premise can be found in the discussion on the application 
limits of local approximation for exchange-correlation interaction of electrons (see, for example, 
[H7] and |38j). On the other hand, an attempt to apply directly an iteration procedure to the 
Kohn-Sham equation with potential U^c in the absence of the Coulomb contribution to the 
effective potential has revealed that iterations failed to converge. 



4.2 Surface energy 

The self-consistent distributions of electron density and potential along with single-particle 
wave-functions obtained in our work allow to calculate the surface energy as a function of Rg. 
However, the formulae for calculating surface energy available in the literature use the finiteness 
of the region occupied by metal. On the one hand, this allows to find the surface energy as 
a difference between total energies of two states of the specimen before and after dividing it 
into two parts. On the other hand, this involves the quasidiscrete spectrum of eigenvalues and 
the necessity of precise evaluation of asymptotic dependence of total energy components on the 
specimen size as it tends to infinity (see, for example, [39]). To calculate the properties of an 
infinite specimen, one should have expressions for the densities of the corresponding components 
of total energy. In the standard formula of the density functional theory for the total energy 
of a finite specimen. 



Etot (A^(r)) = Tg (iV(r)) + J dr6^,[N{r)]N{r) + drdr 



,{ N{v)-N+) {N{v')~N^) 
|r — r'l 



(4.4) 



where integration is implied over the volume, the second and third terms describing the 
exchange-correlation and electrostatic energy of a many-electron system already include den- 
sity definitions of the respective components. As concerns the total kinetic energy Tg of non- 
interacting electron gas, it does not possess such a property, being generally calculated by means 
of the formula obtained in [21] that expresses the sum of single-particle energies till the Fermi 
energy as the integral over the spectrum of the phase of continuous-spectrum wave functions 



with asymptotic form (2.35). Following the formula (50) in [T], we define the kinetic energy 
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density at zero temperature as 



t,(r) = 2 J D{4£|^,(r)|'-f/eff(r)iV(r). 



(4.5) 



e<eF 

Using this expression, one can define the total energy density 

etot(r) = ts(r) + e^,[N{r)]N{r) + ^U{r) {N{r) - iV+(r)) (4.6) 

and, accordingly, obtain the expression for the total energy in the volume V selected out of the 
infinite specimen, 

^tot(V^) = J f/retot(r), (4.7) 

V 

which can be applied for calculating the surface energy. 

The surface energy [3] in the jellium model becomes negative at Rg < 2.5, which means that 
metals with a high density of conduction electrons {2 < Rg < 2.5) could not exist in such model. 
The surface energy values obtained in [3] differ the most drastically from the measured ones for 
multivalent metals, though even the case of monovalent metals of the first group demonstrates 
the increasing worsening of agreement between calculations and experiment with decreasing Rg. 
As it was shown, the main reason is that the jellium model neglects the discrete structure of the 
positive background, which is actually a crystal lattice of ion cores (see, for example, discussion 
in IS])- Nevertheless, it was interesting to find out if the use of self-consistent solution affects 
the Rg dependence of surface energy. 

The surface energy was calculated using two expressions. One of them is based upon the 
strictly thermodynamical definition of surface energy as the component of total energy of a bulk 
homogeneous system, which depends on its surface area S rather than the volume V. Using, 
for instance, the definition (3.82) from |10], one can write 

Etot = 5S + Vetou (4.8) 

where E is the surface energy density. In our case of semi-infinite specimen, the volume of the 
considered part of the system can be modified by receding from the surface to the bulk. Let us 
define the energy per unit area of the system region up to coordinate z by formula 

z 

Etotiz) = [ dzetotiz). (4.9) 



Since the energy density along with all variables defining it becomes constant in the bulk, the 
value of Etot{z) should become a linear function of coordinate at sufficiently large z. Having 
continued this linear dependence back to the boundary of positive background, we obtain a 
finite value independent on the volume, which provides the surface energy by definition. The 
above can be expressed by simple formula 

S= lim [Etot(^) -etot(oo)^]. (4.10) 

2— >00 

It is assumed here that the positive background boundary is at z = 0. 
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Such definition of surface energy is rigorous and, as a proper cliaracteristic of a thermody- 
namical system, independent of the surface-forming process, which is sometimes modeled to 
determine the surface energy. However, it requires an enhanced accuracy in solving differen- 
tial equations and calculating integrals, because the formula includes a difference of two large 



values. This difficulty is eliminated, if Eq. (|4.10|) is rewritten in the form 

S = 



dzttotiz) + dz (ttotiz) - etot(oo)) . (4.11) 



With such a formula for the surface energy, the result proves to be less sensitive to the accuracy 
of calculation, because the integrands tend to zero when receding from the surface. 

In the literature, another definition of the surface energy is generally used, which represents 
it as a difference between the total energy of the semi-infinite specimen and that of a half of the 
infinite one (see, for example, [SI])- Let us represent for convenience, as in work [3], the total 
surface energy density as a sum of kinetic, electrostatic and exchange-correlation components: 

E = Es + Ses + (4.12) 

In the dimensionless form and after integrating over the wave vector parallel to the surface in 



formula (4.5), the components of the surface energy density are written as 



1 

~ 3^ 



dC {HO + HQ]) n{C) - U.C [1] ^(C)} (4.13) 
= ^,J^dC {£.e HO] n{0 - [1] ^^(C)} (4.14) 
= A r ^C{^(C) - 1 ■ OiOMO- (4.15) 



6tt _ _ 

The comparison of Eqs. (4.13)-(4.15) with the formulae (45a)-(45c) in work [34] reveals 
that only the expressions for the kinetic energy contribution ag are different. Using Eqs. (4.5) 
and (4.11) helped to avoid additional complications by eliminating integration of phase (2.37) 
of wave functions (2.35) over energy. 

The results of surface energy calculation are presented in the Table [Tj One can see that 
negative values of the surface energy appear at the values of Rg even slightly larger than those 
reported in work [3]. The main reason for this is a smaller value of correlation energy given 
by formula (2.9). The general conclusion that the jellium model is insufficient to describe the 
surface energy in the domain of existence of real metals is still true. 



5 Calculation of self-consistent potential 
of a Schottky barrier structure 

Up to now, we considered the semi-infinite systems with purely real wave functions. As a 
consequence, such systems do not carry electric current. In this section, a structure is ana- 
lyzed whose single-particle wave functions are limited and show oscillating behavior in both 
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infinitely far regions. Sucli boundary conditions for the Sclirodinger equation solutions admit 
the complex eigenfunctions and system states with the nonzero stationary current. To avoid 
too laborious calculations in actually a model situation, we take a structure such as a metal- 
semiconductor contact with Schottky barrier. Because of a much larger free carrier density in 
the metal the potential barrier is almost totally formed in the semiconductor. This allows one 
to reduce the calculation part of the problem to finding the wave functions and the potential 
in a semi-infinite space and apply to this structure without any essential modification the iter- 



ation algorithm described in Section |2.2| for searching the self-consistent electron density and 
the effective potential barrier. 

We restrict ourselves to the equilibrium situation when no bias voltage is applied to the 
structure and the current is zero. Let the metal occupy the region with z < 0, the degenerate 
n-type semiconductor with the ionized donor concentration A^d located at 2; > 0. For defi- 
niteness, we assume that high density of states at the metal-semiconductor boundary fixes the 
position of the conduction-band bottom at the interface at energy $s reckoned from the Fermi 
level. Since the under-barrier electron density in semiconductor near the interface is small, the 
effective potential coincides with the Coulomb one, thus defining the first boundary condition 
for electrostatic potential. The second boundary condition for solving the Poisson equation in 
the semiconductor region is specified by the requirement for the potential to be bounded at 
infinity, which means at the same time that the electric field turns to zero therein. Reckoning 
the energy from the value of potential in the bulk, we get the boundary conditions for solving 
the Poisson equation in the semiconductor. 



f/(0) = f/0 = $s + £F, 



dU 

dz 







(5.1^ 



Following the common practice, we neglect the penetration of electric field into the metal. 
Thus, the electrostatic potential in the metal is assumed to be constant and equal to its value 
at the interface. 

To describe the single-particle states of continuous spectrum in such one-dimensional inho- 
mogeneous structure, we use the results of the analysis carried out in [2J of the Hilbert space 
related to the Hamiltonian of infinite system. We take as a basis the wave functions of the 
right-hand \Ef^ and the left-hand states 



*fc,k|| = C,^^,^(z)exp (ikiirii) /27r, 

^ik|, (^,r|,) = C,V^(z)exp (ik|,r||) /27r, 

which are uniquely defined by their asymptotic behavior at infinity. 
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(5.2) 
(5.3) 



(5.4) 



(5.5) 



Using two quantum numbers k and q to describe the asymptotic behavior of eigenfunctions 
in the right-hand and left-hand infinite regions is related to the asymmetry of the considered 
barrier structure even in the absence of bias voltage. This asymmetry is caused by different 
positions of conduction-band bottom in the metal and in the semiconductor with respect to 
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the Fermi level. Besides, even a symmetric structure becomes asymmetric when biased, so the 
different quantum numbers are rather a rule than an exception. The relation between k and 
q is found from the condition that both characterize the eigenfunction pertaining to the same 
eigenvalue e of the Hamiltonian, 

Here Uu is the position of conduction-band bottom in the metal reckoned from that in the 
semiconductor, and indices M and S denote the effective masses in metal and in semiconductor, 



respectively. The right-hand equality in Eq. (5.6) defines q as an implicit function of k (and 
vice versa). Normalization constants and equal l/-\/27r if are normalized to 6{q — q') 
and ip^ are normalized to 5{k — k'). The proof of mutual orthogonality of the pair of functions 
^'^'^ at one energy value e, which is necessary to have a full orthonormalized basis in the Hilbert 
space of the problem, was presented in [2]. 

However, in the considered structure, all the results can be expressed in values related to 
the semiconductor electrode of the metal-semiconductor contact, so it is more convenient to 
normalize also to 6{k — k'). In this case, for the normalization constant of the left-hand 
states one should use the expression 

Solving the problems where the self-consistency is unnecessary, one usually normalizes the 
wave functions of continuous spectrum to the unit amplitude of incident wave. It is sufficient 
to calculate the tunneling transparency for the barrier of specified shape. However, to find 
the self-consistent solution taking into account the Coulomb interaction of charge carriers one 
has to know the spatial distribution of electrons, which is impossible without calculating the 
continuous-spectrum eigenfunctions normalized to the delta-function of quantum numbers. 

Let us adduce the relations between the reflectivity r and the transparency t of the barrier, 
which are necessary to carry out the calculations: 

' \tf\' = ^-\rf\\ Jnr\t^\' = l-\r^r (5.8) 



dq/dk I I I I ' dk/dq 



Eqs. (5.9 ) are written for the case of left-hand states normalized to 5{k — k'), which is indicated 
by explicit dependence q{k) in the arguments. 

The electron wave function in semiconductor can be represented as a superposition of two 
real linearly independent solutions (pi and 02, which at 2; — > oo have the following asymptotic 
behavior: 

(pki{z) =sm{kz + 'yk), 
(f>k2{z) = cos{kz + 7fc), 

where 7^ is the phase of the wave function of k-th state. With these functions, using boundary 
conditions and normalization, one can build il'k'i^)^ which, in its turn, is sufficient to define 
^L(^). 

The solutions 0i and 02 are found by means of integration of the Cauchy problem for 



the Schrodinger equation (2.34). With argument z receding from the interface into the bulk, 
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functions 0i increase exponentially in the under-barrier region while 02 decrease. If the po- 
tential barrier is high enough, searching for the exponentially vanishing solution presents some 
calculation difficulties. To overcome them, the following technique was employed. Function 
i)i is found from the Cauchy problem solution with the initial conditions 0fri (0 ) = 1 and 



'/'fci(O) = ^/Uq — specified at the metal-semiconductor interface z = 0. Eq. (2.34) is in- 
tegrated from ;z = to 2;oo into the semiconductor bulk where the effective potential becomes 
constant to a sufficient degree of accuracy. The solution obtained is normalized to the unit 



amplitude by formula similar to Eq. (2.36), and phase 7 in expressions (5.10) is found with 
the help of relation (2.37). It should be noted that choice of initial conditions for 0i is quite 
arbitrary, because the normalization procedure provides the correct asymptotic behavior of so- 
lutions at infinity. To find 02, one also solves the Cauchy problem, but the initial conditions 
are imposed this time at the point z = Zoo'. 

0fc2(^oo) = cos(A;2;oo + 7fc) (5 11) 

0fc2(^oo) = -k(t)kl{Zoo) 

and the integration is carried out from the semiconductor bulk to metal boundary in the 
direction of the (j)k2{z) growth in the barrier region thus improving the stability of procedure. 
The accuracy criterion for the pair of solutions obtained is the constancy of the Wronskian 
W^(0fci,0fe2) at the interval [0,2oo]- 

Further on we assume that the effective masses of electrons in semiconductor and metal are 
equal. Then both the wave function itself and its first derivative are continuous at the boundary 
of the two media. Of course, such premise is rather simulative, but the complicated problem of 
deriving the matching conditions for the envelope wave functions at the metal-semiconductor 
interface is not the subject of our work. A different type of boundary condition at the interface 
may, of course, affect the particular numerical results, but this will hardly require a change in 
the iteration scheme suggested for solving the Kohn-Sham equations for the structures with 
the current states. The reasons for such opinion will be discussed at the end of this section. 
Note also that problems close to those considered here arise as well in the theory of resonant 
tunneling diodes (see, for example, [H]). 



5.1 Some calculation details 



Now then, on the basis of two functions (5.10) we should form the wave function '^^ 



"right-hand" state that describes the tunneling of electrons from the right-hand (semiconductor) 
electrode into the left-hand one, which we still call for convenience "the metal". The asymptotic 
boundary conditions (5.4 ) defining ip^{z) realize the requirement to have a purely outgoing wave 



in the metal and a linear combination of incident and refiected waves in the semiconductor. 
These conditions can be used to find the coefficients in the expansion 



ipf{z) = ak4>ki{z) + Pk4>k2{z) 



(5.12) 



When ipf- is found, in general form can be represented as a linear combination of two 
independent solutions ipf and ipf*. But to calculate the electron density in the semiconductor, 
which is what we only need for the problem to be solved, one has to know ip^iz) only at ^ > 0. 
So to express ipk through 



the left-hand state at z 



bki,2 it is enough to take advantage of the asymptotic form (5.10) of 
00. All the above considering, one easily obtains 



V^g(fc)(^) = t^{k) exp(-i7fc)(i0fci(2;) + 0fc2(^))- 



(5.13) 
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The problem of defining ipfiz) reduces to finding four complex coefficients ak,f3k, t^, t^. The 
four corresponding equations are obtained by equating the asymptotic formulae for function 
ipf' and its derivative at 2; — oo to the asymptotic form of the linear combination required, 
which provides two equations, and also from two matching conditions for the function and its 
derivative at the boundary z = 0. In so doing, we assume that in the metal 



tfexp{-iq{k)z) 



(5.14) 



The calculation is carried out under the condition Uu < 0, which provides that the conduction 
band bottom in the left-hand electrode is lower than that in the right-hand one, so that in- 
equality q{e, k||) 7^ k{e, ky) remains valid in spite of the assumed equality of effective masses in 
the electrodes. At the same time, the condition of effective mass equality ensures that function 
q{k) is independent of ky, which reduces the computational cost. 

With some simple transformations, the sought system of linear algebraic equations takes 
the form 

a0fci(O)+/50fc2(O)-tfc = O, 
a0',i(O)+/50',2(O) + igtfc = O, 
ia + P = 2exp(i7fc), 
ia- P + 2rk exp(-i7fc) = 0. 



(5.15) 



For brevity, we omit everywhere the index R of the right-hand state, the wave-number index k 
of coefficients a , /?, and imply that q = q{k). 



The solution of system (5.15) can be represented as 



rk 



q^k + i$fc 



tk = $fc + rfe$fc, 

ak = -i (exp(i7fe) - rfc exp(-i7fc)) 
Pk = exp(i7fc) + rk exp(-i7fc). 



where for compactness we introduced functions 

= exp(-i7fc) 



and 0'^ means the derivative by z. 



exp(-i7fcj 



^L(0) + i^ 



'fci(O)) 
4i(0)) 



(5.16) 

(5.17) 
(5.18) 
(5.19) 



(5.20) 
(5.21) 



The scheme presented for calculating wave functions describes a modified procedure for 



solving the Schrodinger equation (2.17) for unbounded system as distinct from a semi-bounded 



one (see Section 2.2). The corresponding formula to replace Eq. (2.38) for calculating the 
electron density in a barrier structure is found from the general expression derived in |2]. 
Taking into account that the electron density should be known only at z > 0, and the absence 
of bias voltage, we get 



N{z) = 2 dk rfk|,/ (e(fc,k||)) 



(5.22) 



where f{e) is the Fermi distribution function, and energy eigenvalues e{k, ky) are given by for- 



mula (5.6). The Fermi level position is defined by the neutrality condition in the semiconductor 
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bulk at 2 oo. Using the definitions (5.2)-(5.3), (5.7), and (5.12)-(5.13) along with relations 



(5.8)-(5.9) between the amplitude coefficients of barrier reflection and transparency, we obtain 



from (5.22) an expression for electron density in terms of the two real solutions 4'ki,2{z) found 
numerically. 



,iV(z) 



(27r)3 



dk I rfk,|/ (£(A;,k||)) X 



'Jki 



{z)-i(t)k2{z)\^ . (5.23) 



Here the coefficients r^, a^, and Pk are defined in (5.16)-(5.19). 



The expression (5.23) for electron density provides an explicit representation of the formula 



(2.18) of iteration procedure for the case under consideration. 



5.2 Results and discussion 

Calculations were carried out at the values of parameters typical for the Al/n-GaAs junctions 



[m* = O.OTme, <l>s = 0.9 eV, 



12.5) with ionized donor concentration Nj^ 



10 



18 



cm 



which corresponds to Rs = 0.66. For q, in virtue of the large difference between Fermi energies 
in semiconductor and aluminium, it turned out possible to accept a constant value that was 
chosen equal to the radius of the Fermi surface of electron gas at Rs = 2.07, = 1. Fig. M 




2 - 
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Figure 8: Self-consistent potential for a Schottky barrier structure {Rg = 0.66). Solid line 
represents the potential found with account of exchange-correlation interaction, dashed line is 
the potential in Hartree approximation. The Fermi level position with account of U^c 

is marked. 

represents the coordinate dependence of the self-consistent effective potential calculated with 
account of the exchange-correlation potential of the electrons (curve 1) and in Hartree approx- 
imation (curve 2) near the metal-semiconductor interface. One can see that including Uxc into 
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calculation modifies the potential barrier shape making the drop of potential steeper. This 
difference has a notable impact on the shape of theoretical current-voltage curve of tunneling 
junction, particularly in the case of forward bias when electrons tunnel from semiconductor into 
metal through the effective potential barrier region below the Fermi level (see [TU], |12]-[13])- It 
is demonstrated in the same works that the account of t/xc allows describing more accurately the 
current- volt age curves of real structures. Hence, the representation of the exchange-correlation 
interaction of electrons by means of the local density approximation remains quite accurate for 
the electron states located sufficiently deep under the Fermi level. 

Besides, the use of self-consistent solution for potential in the tunneling current formula 
makes possible to define the parameters and Nj:, of the metal/heavily n-doped GaAs junction 
from the bias voltage dependence of differential resistance |44j-|46Q The Schottky barrier was 
calculated there in the Thomas-Fermi and Thomas-Fermi-Gombas approximations. Though in 
such structures the conditions are well fulfilled for the Schottky barrier to be quasiclassical, 
i. e. k-pL ^ 1 and /cf^tf 3> 1, where L and /tf denote the characteristic scales of spatial 
variations of the self-consistent potential inside and outside the barrier ^44], the question of a 
possible notable change in barrier shape beyond the scope of the Thomas-Fermi approximation 
still remained. The exact calculation showed that the large-scale behavior of the Schottky- 
barrier effective potential is accurately described by the potential found in the Thomas-Fermi 
approximation, and the agreement between them improves with decreasing R^. 




Figure 9: Friedel oscillations of electron density and effective potential in the Schottky barrier 
structure at -Rg = 0.66. Solid line is electron density with account of exchange-correlation 
interaction, dashed line is electron density in the Hartree approximation. Bold solid line is the 
self-consistent effective potential. 

In the semiconductor bulk, the coordinate dependence of self-consistent potential and elec- 
tron density manifests oscillations (Fig. |9]), which are absent in the Thomas- Fermi approxima- 
tion. As a whole, the oscillations caused by Schottky barrier are less pronounced than those in 

^The Fig. 3 missed in the published version of the Ref. [46] is restored in the arXiv preprint. 



34 



the case of an infinitely high potential wall (compare with Fig. |4]), which is apparently related 
to the quasiclassical smoothness of the Schottky barrier. This is particularly true in regard to 
the effective potential oscillations. As can be seen in Fig. |9| the oscillation magnitude of effec- 
tive potential is quite small, of the order of 10~^ and less, while the electron density oscillations 
in the exact solution are, on the contrary, more intensive than those in the Hartree approxima- 
tion. The suppression of the effective potential oscillations compared to those of its Coulomb 
component was already discussed in Section |3.1[ The amplification of density oscillations in 
the exact solution with respect to solution in the Hartree approximation can be explained by 
analyzing the linearized expression for the induced density. 



It is well known that linearization of formula (2.7) in the Hartree approximation at small 
values of potential U leads in the self-consistent Poisson equation (2.12) to the expression for 
the density (charge) increment, which is generally written as 47r(5A^ind = —k^pU. The account 
of the electron density dependence of exchange-correlation potential leads to a similar relation 
4:n6Nin4 = —k^^j.U, the inverse square screening length being 

2 '^TF 
^scr ~ ~^ I 3~1 7 \ / r • (5-24) 

1 + ^aUy,c[n)/dn 

The dimensionless exchange-correlation potential grows in absolute magnitude with decreasing 



electron concentration, so the derivative du^c/dn < 0, and the denominator in formula (5.24) 
tends to zero with growing Rs. The equality of the denominator to zero provides the critical 
value of Rsc- So, at comparable values of the self-consistent potential oscillations in the exact 
solution and that obtained in the Hartree approximation, which is in evidence, the electron 
density oscillations in the exact solution will be amplified increasingly with Rg getting closer 
to -Rsc- Accordingly, the screening region, where the perturbation drops exponentially, will be 
reduced. 

Summarizing the results of analyzing the Friedel-type oscillations caused by the violation 
of homogeneity of degenerate electron gas, we note that the account of exchange-correlation 
interaction of electrons leads to the increase of electron density oscillations and the decrease of 
effective-potential oscillations with respect to results obtained in the Hartree approximation. 

To conclude this section, let us discuss to what extent the results obtained depend on the 
adopted condition of equal effective masses in the two electrodes. Generally speaking, the 
problem of boundary conditions for the envelope function at the interface of two solids differing 
not only by effective masses, but also by the type of charge carriers and their dispersion law, 
is rather complicated. According to the current opinion, such conditions cannot be reduced to 
a universal relation between the values of wave function and its derivative on both sides of the 
interface, but should be derived for each specific case (see, for example, the critical discussion 
related to heteroj unctions in |17], [IE])- 

Staying within the framework of phenomenological description of boundary conditions, one 
can write them in the most general form in terms of the transfer matrix T, the elements of which 
are constrained by the requirement of Hermiticity of the Hamiltonian [18] . If envelope-function 
equations for both electrodes can be taken in the single-band approximation, the boundary 
conditions are described by a 2 x 2 matrix. In the case when the T-matrix can be diagonal 



[19] . the final formulae (5.16)-(5.19) keep their form with the barrier transparency tk and wave 
vector q replaced by renormalized values i = Tutk and q = (T22/Tii)q. 

If matrix T has the nonzero off-diagonal elements, it does not lead to crucial difficulties 
in building the R- and L-solutions of the Schrodinger equation according to their asymptotic 
behavior at z — > ±oo from the particular real solutions. However, the analysis of the most 
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general case in this work would make the statement of fundamental aspects of iteration algo- 
rithm application unreasonably cumbersome. At the same time, the inequality of wave vectors 
q{e) 7^ k{e) typical for the non-symmetrical barrier structures was necessary to include into the 
calculation. 

6 Conclusion 

The reahzed iteration algorithm for solving the self-consistent field equations was found to 
converge in all the cases considered of typical inhomogeneous electron systems with continuous 
spectrum. The difficulties of its application arose only when considering the jclliTim model at 
large values of Rg > 5 that do not occur in most metals and semiconductor structures. It was 
also shown that an explicit account of screening properties of electron gas in the algorithm for 
numerical analysis of infinite Coulomb systems allows to consider a selected fimited volume of 
such system, transferring the boundary conditions from infinity to the surface and assuming 
the energy spectrum of single- electron states to be strictly continuous. Such approach also 
eliminates the difficulties that arise when considering an inhomogeneous infinite system with 
continuous spectrum as a limit of a finite one with quasi-discrete spectrum. 

This work was partially supported by Russian Foundation for Basic Research with grants 
07-02-01481 and 06-02-16955. 



36 



References 

[1] W. Kohn and P. Vashishta, In Theory of the inhomogeneous electron gas, ed. by S. 
Lundqvist and N.H. March, Plenum Press, N-Y. (1983), Ch. 2 



[2] 
[3] 
[4] 
[5] 

[6] 
[7] 
[8 
[9 
[lO: 

[11 
[12 

[13 

[14 

[15 

[16 

[17 
[18 
[19 
[20 
[21 
[22 
[23 
[24 
[25 



A.Ya. Shul'man, J. Phys.: Conf. Ser. 35, 163 (2006) 
N. D. Lang and W. Kohn, Phys. Rev. B 1, 4555 (1970) 

M. Manninen, R. Nieminen, P. Hautojarvi, J. Arponen, Phys. Rev. B 12, 4012 (1975) 

A. Liebsch, Electronic Excitations at Metal Surfaces, Plenum Press, New York (1997), Sec. 
2.3.2 

R.M. Nieminen, J. Phys. F, 7, 375 (1977) 



A.Ya. Shul'man, D.V. Posvyanskii, arXiv: cond-mat/0209335 



D. Pines, Elementary Excitations in Solids, Benjamin, N-Y, 1963 

P. Gombas, Theorie des Atoms und Ihre Anwendungen, Springer Verlag, Wien (1949) 

A.Ya. Shul'man, V.V. Zaitsev, Sol. State Comm., 18, 1623 (1976) 

S.K. Godunov, V.S. Ryaben'kii, Difference schemes, Nauka, Moscow (1977) (in russian) 

O.V. Konstantinov, A.Ya. Shik, ZhETF, 58, 1662 (1970) 

J.F. Appelbaum and G.A. Baraff, Phys. Rev. Lett 26, 1432 (1971) 

R.G. Newton, Scattering Theory of Waves and Particles, McGraw-Hill, New York, 1966 

G.A. Baraff and J.F. Appelbaum, Phys. Rev. B 5, 475 (1972) 

R.P. Fedorenko, Introduction in Computational Physics, Izd. MFTI, Moscow (1994), §15 
(in russian) 

T. Ando, A.B. Fowler, F. Stern, Rev. Mod. Phys. 54, 437 (1982) 
J. Suiie, P. Olivo, B. Ricco, J. Appl. Phys. 70, 337 (1991) 

G. A. Mead, Phys. Rev. Lett. 6, 545 (1961) 

M. Stengel and N.A. Spaldin, Nature 443, 679 (2006) 

H. Y. Ku and F.G. Ullman, J. Appl. Phys. 35, 265 (1964) 

G. T. Black and J.J. Welser, IEEE Trans. ED 46, 776 (1999) 
J. Bardeen, Phys. Rev. 49, 653 (1936) 

H. B. Huntington, Phys. Rev. 81, 1035 (1951) 

L.J. Sinnamon, R.M. Bowman J.M. Gregg, Appl. Phys. Lett. 78, 1724 (2001) 



37 



[26] H. Ehrenreich and H.R. Philipp, Phys. Rev. 128, 1622 (1962) 

[27] W.S. Choi, S.S.A. Seo, K.W. Kim et all, Phys. Rev. B 74, 205117 (2006) 

[28] K.M. Rabe, Nature Nanotechnology 1, 171 (2006) 

[29] W. Kohn and C. Majumdar, Phys. Rev. 138, A1617 (1965) 

[30] D. C. Tsui, Phys. Rev. B 8, 2657 (1973) 

[31] N.D. Lang and W.Kohn, Phys. Rev. B 3, 1251 (1971) 

[32] J. P. Perdew and Y. Wang, Phys. Rev. B 38, 12228 (1988) 

[33] H.F. Budd and J. Vannimenus, Phys. Rev. Lett. 31, 1281, 1430 (Erratum) (1973) 

[34] N.D. Lang, In Theory of the inhomogeneous electron gas, ed. by S. Lundqvist and N. H. 
March, Plenum Press, N-Y. (1983), Ch. 5 

[35] B. Sieroszynska-Wojas, FTT 10, 693 (1968) 

[36] R. Monnier and J. P. Perdew, Phys. Rev. B 17, 2595 (1978) 

[37] R. O. Jones and O. Gunnarsson, Rev. Mod. Phys. 61, 689 (1989) 

[38] W. Kohn, Rev. Mod. Phys., 71, 1253 (1999) 

[39] G. Paasch and M. Hietschold, Die Oberfldche der festen Korper. In: Ergebnisse in der 
Elektronentheorie der Metalle , eds P. Ziesche and G. Lehmann, Akademie-Verlag, Berlin 
(1983), Ch. 10 

[40] M.A. Leontovich, Introduction in Thermodynamics, Nauka, Moscow (1983), Ch. 3 (in 
russian) 

[41] D. K. Ferry and S. M. Goodnick Transport in Nano structures, Cambridge Uni Press, 
Cambridge, UK (1997) 

[42] I.N. Kotel'nikov, A.Ya. Shul'man, Proc. 19th Int. Conf. on Phys. Semicond., Warsaw 
(1988), Vol. 1, p. 681 

[43] A.Ya. Shul'man, I.N. Kotel'nikov, N.A. Varvanin et al, Pis'ma v ZhETF, 73, 643 (2001) 

[44] I.N. Kotel'nikov, I.L. Beinikhes, A.Ya. Shul'man, FTT 27, 401 (1985) 

[45] I.N. Kotel'nikov, D.K. Chepikov, E.G. Chirkova et al, FTP 21, 1854 (1987) 

[46] E.M.Dizhur, A.Ya. Shul'man, I.N. Kotel'nikov, A.N.Voronovsky, phys. stat. sol. (b) 223, 
129 (2001); |arXiv:cond-mat/0010200| 

[47] E.E. Takhtamirov, V.A. Volkov, ZhETF 116, 1843 (1999). 

[48] I.V. Tokatly, A.G. Tsibizov, A.A. Gorbatsevich, Phys. Rev. B 65, 165328 (2002) 
[49] T. Ando, S. Wakahara, H. Akera, Phys. Rev. B 40, 11609 (1989) 



38 



